A self-organized critical model of rearranging hydrogen-bonded network in ice 
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A dynamical four- vertex model of rearranging hydrogen-bonded network in ice, with a dynamics 
similar to the Grotthuss mechanism, is proposed. The model qualitatively explains the unusually 
high proton-conductivity in ice. It also serves as an interesting example of self-organized criticality 
with a non-conservative dynamics. The model is solved on a linear chain and its steady state is 
determined. Using numerical simulations, the model is further studied on a square lattice, and it 
is observed that it relaxes by three different avalanches, each of which follows a simple finite-size 
scaling. 
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An isolated water molecule can be said to be well un- 
derstood because its properties, as observed in experi- 
ments, can be described by the first principles quantum 
mechanical calculations. However, when a large number 
of these molecules form liquid- water or ice, our ability 
to predict their properties becomes limited In fact, 
water has a large number of anomalous properties Q, 
amongst which the unusually high mobility of H + has 
attracted a lot of attention. The observed value of the 
mobility in bulk ice 0] or water-filled narrow pores [H 
are comparable to the electronic mobilities in some semi- 
conductors. Conduction of H + is also fundamental to a 
myriad of processes ranging from ATP synthesis to elec- 
trical power generation in hydrogen fuel cells. 

The generally accepted idea is that the hydrogen- 
bonded network of water molecules provides conduction 
pathways to H + , which migrates via Grotthuss mecha- 
nism [5[. In this, an extra H + in the network appears 
as a hydronium ion (H30 + ), which translocates by hop- 
ping of H + across the hydrogen bonds (Fig. [T]). This 
movement is hindered by the broken hydrogen bonds in 
the network (now on referred to as defects), which them- 
selves migrate by rotation of the water molecules. The 
rotation of the molecules is slower than the hopping of 
H + and is therefore the rate limiting step However, 
a comparison of the mobility of H + , determined experi- 
mentally, with that of the defects, estimated assuming a 
simple diffusion, shows that the former is much higher at 
low temperatures 0] . We show that this problem can be 
qualitatively resolved if one considers a long-range corre- 
lation in the hydrogen-bonded network. While presence 
of this correlation has been experimentally observed and 
its importance envisioned 0, a clear understanding is 
still missing. In this letter, we propose a model of re- 
arranging hydrogen-bonded network, with a Grotthuss- 
type dynamics, which exhibits a long-range correlation in 
its steady state. This correlation in the network increases 
the mobility of the defects, and a long distance transport 
of these occurs with a non-zero probability. We solve 
this model on a linear chain and determine its steady 



I 



v. 



D I 



I 

HP* V 



L I 



H ^ 



I 



r 



'o 
I 



cr 
I 



FIG. 1. (Color online) A schematic of H + transport by the 
Grotthuss mechanism. The hydronium ion and the two types 
of broken bonds, D and L are shown by red letters. The 
hopping of H + and the rotation of molecules around an O-H 
axis are indicated by arrows. 



state. Further, using numerical simulations on a square 
lattice, we show that, in the steady state, the probability 
distribution of long distance transport of the defects can 
be described by a linear combination of simple scaling 
forms. 

The model is interesting per se, as it provides a 
theoretical example of Self-Organized Criticality (SOC) 
with non-conservative dynamics. The existence of non- 
conservative models of SOC has been long debated. The 
much studied models, like the forest fire model @, the 
Olami-Feder-Christensen model in the non-conservative 
limit 0] and their variants are known to be non-critical 
[loj . The non-conservative sandpile model studied by 
Prucssner and Jensen ll| is critical in the infinite vol- 
ume limit, but lacks a well defined finite-size scaling, un- 
less a control parameter is fine-tuned ■ The model we 
propose does not have any conserved quantities. When 
slowly driven, it exhibits a scale invariance, following a 
finite-size scaling, without any apparent fine tuning of 
parameters. As the model shares some features with the 
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FIG. 2. An illustration of the relaxation, started by rotating 
the molecule at the top-left site of a 2 x 2 lattice. The covalent 
O-H bonds are denoted by solid lines, and the rotation of 
molecules are indicated by arrows. The four orientations of a 
molecule are denoted by 1, 2, 3 and 4. 



well known ice model [12j of equilibrium statistical me- 
chanics and also with the sandpile model of SOC [ill, we 
call it the icepile model. 

We define the model on a square lattice with open 
boundaries (Fig. [2]). Every site on the lattice is occupied 
by an oxygen atom. Two hydrogen atoms are connected 
to each site by covalent bonds and they are always at 
right angles to each other. This is a special case of the 
standard ice model, with only four out of the six orienta- 
tions allowed, i.e., the orientations where two hydrogens 
are diametrically opposite to each other, are not allowed. 

Any orientational arrangement of the molecules may 
have empty edges or edges occupied by two hydrogen 
atoms. These are known as L- and D-defects [Bj], respec- 
tively. Hydrogens in the neighboring molecules interact 
via repulsive electrostatic interaction, and the D-defects 
cost energy. Let the energy be U for a D-defect, and zero 
for the singly- or un-occupied edges. Further, to specify 
the orientations of the water molecules, we introduce a 
set of variables Tij such that Ti j — 1, when a hydro- 
gen atom from the ith site is along the i-j edge, and 
Ti.j = 0, otherwise. The energy of a configuration is then 
described by an effective Hamiltonian 
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where the angular brackets indicate that i and j are near- 
est neighbors, and Tjj = for sites outside the boundary. 
The equilibrium properties corresponding to this Hamil- 
tonian are easy to determine, if we represent each orien- 
tation of a water molecule by a pair of Ising-spins, taking 
up or down and left or right orientations. This breaks the 
system into mutually uncoupled horizontal and vertical 
linear chains of Ising-spins, which is exactly solvable. 

As mentioned before, the translocation of the defects 
determines the effective mobility of H + , and this is what 
we study in the non-equilibrium steady state of our 
model. We assume the configurations with D-defects as 



unstable. In an ice crystal, the defects are formed mostly 
by thermal excitations and their diffusion takes place 
by the reorientation of water molecules We gener- 
ate the defects by forcibly reorienting a randomly chosen 
molecule by 90° clockwise rotation, calling it an exter- 
nal drive. If this creates a D-defect, we use the following 
odd-even update rule to relax the configuration (see Fig. 
[2]): We divide the lattice into two sub- lattices according 
to the parity of x + y, where x and y are the Cartesian co- 
ordinates of a site. We pick one sub-lattice and select all 
sites on it contributing to a D-defect and simultaneously 
rotate the corresponding water molecules by 90° clock- 
wise. This may create or annihilate D-defccts and/or 
move them by one lattice unit. We do the same for 
the other sub-lattice and repeat until the system reaches 
a stable configuration. At this stage we again disturb 
the system by forcibly rotating another randomly chosen 
molecule and repeat. 

We refer to a complete relaxation process, starting 
from the driving to reaching a stable configuration, as 
an avalanche. The number of D-defects is not conserved 
in an avalanche. The reorientation of the molecules in a 
relaxation step may create more D-defects which move on 
the lattice until they meet an empty edge or the lattice 
boundary, where they annihilate. Unlike the conserved 
sandpile models, the activity during a relaxation is not 
coupled to any conserved quantities. 

There are a large number of stable configurations, all 
of which appear with a non-zero probability in the steady 
state. This can be verified from the fact that the config- 
uration with all water molecules in the same orientation 
can be reached from any stable configuration, and vice 
versa. In a stable configuration, each row and column in 
the lattice has at most one empty edge whose positions 
uniquely determine the configuration. For each row with 
L sites, there are L + 1 ways of placing the empty edge 
(we consider the configuration with no empty edges as 
the one with empty edges outside the lattice). Then for 
an L x M lattice, there are 



(L + 1) M (M+l) 1 



(2) 



ways of placing the empty edges, thus, so many stable 
configurations. 

We first study an L x 1 chain and show that all the 
2 L (L + 1) stable configurations are equally probable in 
the steady state. Let |C) be a stable configuration. De- 
fine an operator a; by the relation ai\C) = \C), where 
\C) is the configuration reached by driving at site i of the 
chain in configuration \C), and relaxing. The D-defect 
created by the driving, moves toward the empty edge and 
vanishes by annihilating with it. It is always the molecule 
at the site next to the empty edge, which rotates last. For 
the present case, let that site be j. We define an operator 
R, corresponding to the reflection of the system across a 
line parallel to the chain, which maps any allowed config- 
uration to another unique allowed configuration. Then it 
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can be shown that there exists an avalanche correspond- 
ing to the equation ajR\C) = R\C). Similar argument 
applies for all the driving operators, implying that they 
all have an inverse. Thus the Markov evolution satisfies 
the pairwise balance and hence the steady state is 
equally probable. 

Let s be the total number of 90° rotations of the water 
molecules in an avalanche, and P (s, L) be the probability 
of such events on a linear chain of L sites. As noted, 
the D-defect generated by the driving moves toward the 
empty edge, and vanishes by annihilating with it. All the 
sites in its path rotate at least once. If the driving site is 
on the right of the empty edge, then, before relaxation, 
the sites in between are either in state 3 or 4 (Fig. [2]). 
For the configurations where all these sites are in state 
3, the corresponding water molecules rotate just once. 
If any of these sites are in state 4, then, each of them 
individually increases the total number of rotations by 
2. Thus, effectively, the sites in state 3 and 4 contribute 
,s = 1 and s = 3, respectively. A similar argument applies 
for the driving sites on the left of the empty edge. From 
this it is easy to show that the probability 
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where [^J is the largest integer less than or equal to x 
and Sij is the Kronecker delta, which corresponds to the 
fact that half of the times the driving does not create 
any D-dcfcct. For a large s, the expression reduces to a 
simple scaling form 



P(s,L)~L- l g(s/L), 



(4) 



where g{x) = 2" 1 (1 - x/2). 

Let us now consider an L x L lattice. Our numerical 
studies show that, in the steady state, the empty edges 
are more likely to be found near the center of the lattice. 
This implies that the stable configurations in the steady 
state are not equally probable. 

On this lattice there are three types of avalanches. A 
schematic of these avalanches are shown in Fig. [3] In 
type I, the activity is confined along a single row or a 
column. In the other two types, the activity spreads on a 
two-dimensional area. In one of these, the activity moves 
like a single wave crossing the sites in its path only once, 
and the water molecules rotate less than five times. This 
is the type II avalanche. In the other, the waves move 
back and forth on the lattice, until the system reaches 
a stable configuration. This is the type III avalanche. 
The number of rotations increases for the sites which are 
away from the boundary of the affected region, and has 
a maximum value of the order L. Then the cutoff in the 
total number of rotations, s, in an avalanche, scales as L, 
L 2 and L 3 for the three types of avalanches, respectively. 
In such a case, the probability distribution is expected to 
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FIG. 3. (Color online) A schematic of the affected regions in 
the three types of avalanches on a 20 x 20 lattice. Each colored 
box is a unit cell of the lattice. Different colors represent, 
different number of rotations per site. 
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FIG. 4. (Color online) The scaled size distribution plots of 
the three avalanches. The exponents for the type I, II and III 
avalanches are a = 1, 2, 3, and U) = 9/8, 11/4, 3, respectively. 
For a better view, the graphs for type II and III have been 
shifted by 10 5 and 10 12 units, on the x-axis. 



have a scaling form 
1 



P(s,L) 
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with the scaling functions /, g and h, which vanish as 
their argument approaches 0(1). A similar linear com- 
bination of simple scaling forms was found in a few other 
models [TBI, EH . 

We have analyzed our Monte-Carlo simulation results 
for five different values of L, each averaged over 10 6 tri- 
als. The best data collapse of the numerical results is 
obtained for the values a ~ 9/8, (3 ~ 11/4 and v ~ 3 
(Fig. [4]). This implies that the probability of type I and 
type II avalanches decreases as L -1 / 8 and L~ 3 / 4 , for large 
L. This is consistent with the numerical results in Fig. [5] 
Thus, in the thermodynamic limit, there will be mostly 
type III avalanches. 

The geometry of the lattice is important in deciding the 
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FIG. 5. (Color online) Probability of different types of 
avalanches. The fitting functions are w (x) = 1.04/x 1 ^ 8 , 
x( x ) = 2.85/x 3//4 and n{x) = 1 — rj (x) — x{x), for type 
I, II and III, respectively. 



criticality of the steady state. To see this, we consider a 
straight-forward generalization of our model on to other 
lattices. First we consider a three-coordinated Caley tree. 
Using an argument similar to the one-dimensional chain, 
one can see that all the stable configurations are equally 
probable in the steady state. Then using the tree struc- 
ture of the lattice, it is easy to show that the probability 
of an avalanche in which r distinct water molecules ro- 
tate, has a form 2 r exp (— Ax/ logy), where v is the num- 
ber of vertices on the Caley tree and A is a numerical con- 
stant. Then the steady state is not critical. On any other 
three-coordinated lattice, e.g., hexagonal lattice, there 
will be lesser freedom for the water molecules. However, 
on lattices like square or Kagome lattice, where each site 
is connected to four neighbors, the probability distribu- 
tion has a power-law tail, implying that the steady state 
is critical. On a lattice with coordination number greater 
than four, e.g., cubic lattice, the water molecules become 
stable more easily, and the probability of getting large 
avalanches decreases exponentially. 

In an ice crystal, the water molecules are arranged on 
a four-coordinated diamond-like lattice. Such a lattice 
is more like the Kagome lattice, and the steady state of 
the icepile model is expected to be critical. Although 
the dynamics of the model in its present form is not very 
realistic (a more careful modeling of the intra-molecular 
interaction and the dynamics is required), it does show 
an important qualitative property of ice - there is a long- 
range correlation in the hydrogen-bonded network Q. 

In the Grotthuss mechanism, because the molecular re- 
orientation is slower than the H + hopping, the effective 
mobility of H + is limited by that of the defects. How- 
ever, when the defect-mobility was estimated^ Q using 
experimental value of water rotation time 



less compared to the experimental value of H + mobility, 
at low temperatures. This discrepancy can be resolved 
qualitatively, staying within the Grotthuss formalism, if 
one considers the long-range correlation in the hydrogen- 
bonded network. The correlated walk on the network in- 
creases the mobility of the defects significantly, and this 
in turn increases the mobility of H + . 

In summary, we have proposed a dynamical model of 
the hydrogen-bonded network in ice, the steady state of 
which, in the slow driving limit, has a long-range correla- 
tion. Ignoring this correlation in the study of proton con- 
duction in real ice, can be misleading. The model is also 
an example of SOC with non-conservative dynamics. An- 
other important aspect of the model is that its dynamics 
is non-abelian, i.e., the final configuration depends on the 
order in which the system is driven. All the known undi- 
rected models of SOC with exact solutions are abelian 
1 3| . There are only a handful of analytical results avail- 
able for the non-abelian models [HI, EH , but those are 
all in the special limit, where the dynamics essentially 
becomes abelian. In fact, even the one-dimensional non- 
abelian models also are highly non-trivial. We solved the 
icepile model in one dimension, and determined its steady 
state. We also calculated the exact number of recurrent 
configurations on a square lattice. A similar calculation 
is possible on a Kagome lattice, and it can be shown that 
the number of recurrent configurations 



n Kagome = V/x4 2L {2L + l) 



(6) 



on a 2L x 2L Kagome lattice lattice with open boundaries. 
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